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LAGRANGIAN  SIMULATION  OF  TAYLOR-COUETTE  FLOW 


I.  Introduction 

We  report  here  on  the  development  of  the  hydrodynamics  code  SPLASH 
which  is  designed  for  the  Lagrangian  simulation  of  transient  rotational  flow 
phenomena.  The  code  solves  the  incompressible,  inviscid  fluid  equations  in 
an  axisymmetric,  cylindrical  coordinate  system.  This  is  a  2  1/2  dimensional 
model  with  two  spatial  coordinates  (r,z)  and  three  velocity  coordinates 
(u,w,v).  All  variables  are  independent  of  the  angle  0. 

The  equations  of  motion  are  f inite-dif f erenced  on  a  general-connectivity 
triangular  mesh.  A  triangular  mesh  is  the  natural  choice  for  flows  in 
complicated  geometries  or  flows  with  free  surfaces  or  interfaces.  The 
general-connectivity  mesh  allows  local  mesh  restructuring  whenever  the  grid 
distorts  sufficiently  to  affect  numerical  accuracy  and  convergence.  The 
SPLASH  code  is  a  direct  extension  of  the  hydrodynamics  code  SPLISH^  which 
solves  the  incompressible,  inviscid  hydrodynamic  equations  on  a  triangular 
mesh  in  Cartesian  geometry. 

This  model  is  applicable  to  a  host  of  problems  such  as  rotating  columns 

of  fluid  including  imploding  liner  systems  with  axially  displaced  annular  or 
2 

end  plate  pistons  ,  axisymmetric  jets,  laser  ablation  of  spherical  shells 

and  droplet  combustion.  Here  we  apply  the  model  to  the  study  of  Couette  flow 
3 

and  Taylor  vortex  formation  between  two  rotating  coaxial  cylinders.  This 

problem  was  chosen  as  a  test  case  because  the  linear  theory  is  straight- 

4  5 

forward  and  well-developed  ’  and  there  is  a  myriad  of  experimental 
6  7  8 

results  ’  ’  available  for  comparison.  We  have  been  unable  to  find  any 
numerical  work  in  the  literature  which  models  the  time-evolution  of  rota¬ 
tional  flows  to  serve  as  a  comparison. 
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In  Che  next  section  we  discuss  several  aspects  of  the  triangular 
gridding  techniques.  The  equations  of  motion  are  developed  in  Section  Ill 
with  the  finite-difference  algorithms  presented  in  Section  IV.  In  Section  V 
we  discuss  the  Couette  flow  problem  and  the  numerical  results  are  presented 
in  Section  VI.  The  summary  and  conclusions  make  up  Section  VII. 


II.  Triangular  Mesh 

The  set  of  equations  governing  incompressible,  inviscid  flow  in  a 
cylindrical  coordinate  system  will  be  approximated  by  finite  differences 
on  a  triangular  mesh.  The  variables  in  these  equations  will  be  represented 
as  triangle  or  vertex  quantities  on  this  mesh.  This  differencing  procedure 
is  somewhat  complex.  We  will  illustrate  it  by  discussing  some  important 
basic  concepts. 

The  basic  computational  cell  is  the  shaded  region  shown  in  Fig.  1.  It 

is  formed  by  joining  the  side  bisectors  of  the  triangles  surrounding  the 

general  vertex.  Each  triangle  surrounding  a  central  vertex  contributes 

1/3  of  its  area  to  the  area  of  the  basic  cell. 

We  now  illustrate  how  a  gradient  is  represented  in  this  model.  If 

vertex  quantities  are  linear  functions  of  position,  then,  given  the 

function  g  (defined  on  vertex  m) ,  the  function  g  at  any  other  point,  n, 
m 

say,  can  be  written  without  approximation,  as 


8n  = 


m 


+  R 
—  n 


(II-l) 


Here  R  is  the  vector  from  the  location  of  g  to  the  chosen  point. 

—  n  m 

Now  consider  the  triangle  j  defined  by  two  side  vectors  ,  _S 
with  the  vertex-defined  quantities  g,  g^  and  The  index  j  indicates 

triangle  quantities,  the  index  i  vertex  quantities  (see  Fig.  2).  Following 
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Fig.  1  —  Detail  of  the  triangular  grid  elements.  A  triangle  is  made  up  of  three 
directed  line  segments,  and  a  vertex  represents  the  shaded  region  defined  by 
the  center  of  mass  of  each  surrounding  triangle  and  the  midpoint  of  each  side. 


* 


Ref.  9,  the  gradient  of  g,  uniquely  defined  on  the  triangle  and  constant 
throughout  it  is, 


Vg.  -  (grg)  s±+1 


(gi+l-g)  i. 


■  I  n 

i=*l 


Ik. 

J 


(II-2) 


where  S_  is  the  side  vector  rotated  clockwise  by  tt/2  radians,  and  is  short¬ 
hand  for  a  cross  product  with  a  basis  vector.  Here  n  is  a  unit  vector  normal 
to  the  computational  plane  and  is  the  area  of  the  triangle. 

2Aj  =  (jr  i  x  (£.“  JL  i+^)  •  n.  _r^  is  the  coordinate  location  of  the 
vertex  i.  The  conclusion  is  that  gradients  may  be  naturally  represented  on 
triangles,  and  easily  calculated  on  them  in  the  linear  approximation. 

The  integral  operator  consists  of  a  piecewise  triangle  summation 
about  a  basic  cell  giving  rise  to  a  vertex  integral  which  is  likewise  exact 
for  the  linear  approximation.  Expressions  for  the  divergence,  curl  and 
Laplacian  are  presented  in  detail  in  Ref.  1.  The  finite  differencing  of 
the  nonlinear  diffusion  equation  (Eq.  A-5)  is  discussed  in  Ref.  10.  Although 
the  basic  difference  and  integral  operators  are  linear,  the  resulting 
weighting  to  a  central  vertex  yields  second  order  accurate  approximations 
in  the  same  way  that  central  differences  are  second  order  accurate  for  one- 
dimension. 

Although  the  control  volume  approach  assures  that  the  equations  are 
solved  conservatively,  large  numerical  errors  may  arise  in  a  Lagrangian  code 
due  to  severe  grid  distortion.  If  portions  of  the  grid  become  stretched, 
gradients  will  be  calculated  which  involve  vertices  far  removed  from  one 
another.  The  convergence  of  the  iterations  would  be  slow  and  truncation 
errors  would  build  rapidly  as  the  triangle  sides  lengthen. 


This  difficulty  is  avoided  by  forcing  the  mesh  to  restructure.  Mesh 
restructuring  may  involve  interchanging  the  diagonals  of  a  quadrilateral 
formed  by  two  adjacent  triangles  or  adding/deleting  a  vertex  on  a  triangle 
side  or  in  the  interior  of  a  triangle.  Restructuring  is  performed  to  pre¬ 
serve  the  diagonal  dominance  of  the  Poisson  equation  as  a  specific  condition 
on  the  representational  accuracy . 

We  employ  the  same  basic  restructuring  algorithms  as  developed  by 
Fritts^  for  Cartesian  geometry.  By  conserving  the  linear  momentum,  circula¬ 
tion,  divergence  and  the  r  and  z  components  of  the  angular  momentum  on  a 
quadrilateral  (or  triangle)  a  unique  and  reversible  solution  for  the  new  r 
and  z  components  of  the  triangle  velocities  is  obtained.  Conserving  the  0 
component  of  the  angular  momentum  and  the  generation  of  circulation  gives 
rise  to  a  unique  solution  for  the  new  angular  momentum.  These  reconnection 
algorithms  will  be  discussed  in  more  detail  in  a  future  report.  In  the  next 
section  we  discuss  the  equations  of  motion. 


III.  Equations  of  Motion 

The  hydrodynamic  equations  of  motion  governing  an  incompressible, 
inviscid  fluid  in  a  cylindrical  coordinate  system  (r,  9,  z  with  corres¬ 
ponding  velocity  components  u,  v,  w)  are 

3u  ,  3u  ,  3u  v2  1  3p 

3t  3r  3z  r  p  3r 


(III-l) 


3v  ,  3v  .  3v  .  vu  . 
— -  +  U  —  +  w  —  +  —  =  0 
3t  3r  3z  r 


(III-2) 


3w  ,  3w  3w  13p 

—  +  u  —  +  w  —  =■ - , 

3t  3r  3z  p  3z 


along  with  the  incompressibility  condition 


(II 1—3 ) 
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(III-4) 


i  a  ,  .  .  a  A 

7  7F  (ru)  +  a7  w  3  0 


Equation  III-2  is  just  the  conservation  of  angular  momentum  per  unit 
mass.  Defining  the  angular  momentum  per  unit  mass  L  =  vr ,  Equation  III-2  is 


3L  ,  m  dL 
— +u-  VL=— =0 


(III-5) 


We  define  a  pseudo-Cartesian  vector  space  with  a  gradient  operator 

,  (III-6) 


rr  *  -  *  3  ,  «  3 

V '  =  e  t—  +  e  t 
r  3 r  z  3z 


a  vector  velocity 


u'  =  e  u  +  e  w 
—  r  z 


and  a  position  vector 


( II 1—7 ) 


r'  =  e  r  +  e  z 


(III-8) 


Equation  III-l  and  III-3  can  then  be  combined  to  give 


,  d  V  '  , 

(dc>  H. 


1  ~ 
—  V  p  +  —  e 
p  r  r 


where 


(III-9) 


<£>’  +  V 


Rewriting  the  r-component  of  Eq.  III-9  we  have 

L2  3£ 
or  ~p7T  =  ~  ^ 

This  is  precisely  the  equation  of  motion  that  would  be  obtained  from  a  one 
dimensional  Lagrangian 

1  1  L" 


V  1.2  1  L 

=■  2  pr  “  2  *7*  ~ 

with  a  potential 


(III-10) 


7 


K  fTP 


(III-ll) 


V 


p 


L2 

pr2 


The  equations  of  motion  in  our  pseudo-Cartesian  coordinate  system  become 

(£)’l-° 

and  the  angular  momentum  of  a  fluid  element,  per  unit  mass,  remains  constant 
as  we  follow  it  with  its  motion;  and 


<-{cy-  =  ■  0  J'p  "  I  l2v'  r*  •  (HI-12) 

The  radial  and  axial  motion  takes  place  as  if  v  were  absent  and, 
instead,  a  centrifugal  force,  =  -  py  L2V'  ,  were  acting  in  the  radial 

direction.  Note  that  the  angular  momentum  is  constant  in  the  region  over 
which  the  gradient  is  applied.  It  is  therefore  a  triangle  quantity  —  the 
basic  fluid  element  in  this  code  —  and  each  triangle  resides  in  a  i/r2 
potential.  The  incompressibility  condition  takes  the  form 

V'  •  (ru’)  =  0  .  (III-13) 

Equations  III-ll,  12,  13  now  have  the  same  formulation  as  the 
Cartesian  version  of  the  code  (SPLISH) .  The  finite-difference  algorithms 
and  the  restructuring  algorithms  developed  for  SPLISH  can  then  be  applied 
to  cylindrical  geometry  with  modifications  to  include  conservation  of 
vorticity  generation  and  conservation  of  angular  momentum.  The  pseudo- 
Cartesian  coordinate  system  will  be  utilized  throughout  the  remainder  of 
this  paper.  For  that  reason  the  primed  notation  will  be  dropped  in  what 
follows.  In  the  next  section  we  discuss  the  finite-difference  algorithms. 
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IV.  Finite-Difference  Algorithms 


The  computer  code  SPLASH  is  a  2  1/2-dimensional  Lagrangian  fluid 
dynamics  code  for  incompressible  fluids  in  cylindrical  geometry.  This  code 
is  a  direct  extension  of  the  philosophy  and  numerical  techniques  developed  for 
SPLISH,  the  Cartesian  version  of  the  code.  As  such,  most  remarks  in  this 
section  apply  equally  as  well  to  SPLISH. 

The  basic  equations  are. 


(IV-1) 


or  explicitly  if  velocities  are  considered  to  be  triangle-centered.  This 
placement  of  velocities  as  cell  quantities  and  pressures  at  vertices  is 
apparently  unique  to  SPLASH  and  SPLISH  and  is  the  direct  opposite  of  the 
usual  placement.  In  what  follows  the  subscript  i  will  denote  a  vertex- 
centered  quantity  and  j  a  triangle-centered  quantity.  In  both  codes  the 
integration  of  velocities  uses  a  split  step  algorithm  whereby  the  velocities 
are  advanced  one  half  timestep  (Eq.  IV-4) ,  the  grid  is  advanced  a  full 


(IV-6) 


R^  *  +  5tu  * 


0  j  =  R[(R°)  ,(rJ)]  •  0  l1 


(IV-7) 


^2-^(7p)rH(7K 


(IV-8) 


The  vertex  velocity  U_ ^  appearing  in  Eq.  IV-5  is  obtained  from  the  area- 


weighted  .  from  the  previous  iteration, 

IunA. 

■  ~J  J 
U  n  =  1-—--- 
-i  EA. 

J 


(IV-9) 


The  advantage  of  using  triangle  centered  velocities  is  the  ease  in  concep¬ 
tualizing  and  expressing  conservation  laws.  Because  of  the  paucity  of 
experience  in  formulating  algorithms  over  a  general  triangular  grid,  we 
employed  a  control  volume  approach,  which  uses  an  integral  formulation  to 
derive  the  difference  algorithms.  Equation  IV-7  is  the  first  manifestation 
of  this  approach.  It  reflects  numerically  the  fact  that  the  triangle 
velocities  must  rotate  and  stretch  as  the  grid  rotates  and  stretches. 

The  transformation  R  is  derived  by  considering  the  circulation  about  each 
vertex.  The  boundaries  of  a  vertex  cell  are  defined  by  the  triangle  side 
bisectors  as  noted  in  Section  II. 

The  vertex  of  Fig.  1  is  constructed  by  summing  over  all  the  surrounding 
triangles.  Therefore  the  area  of  a  vertex  cell  may  be  defined  as 


A  =Ita-  ,  (IV-10) 

c  .  3  j 
J 

where  the  sum  extends  over  all  adjacent  triangles.  With  this  definition  the 
vertex  velocity  becomes 
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c 


Since  the  triangle  velocities  are  constant  over  the  triangle,  the  circula¬ 
tion  taken  about  the  boundary  of  the  vertex  cell  is  straightforward.  Circu¬ 
lation  is  conserved  about  each  of  the  three  triangle  vertices  by  the  trans¬ 
formation  R  of  Eq.  IV-7.  This  transformation  ensures  that  the  vorticity 
integral  calculated  about  any  interior  vertex  is  invariant  under  the 
advancement  of  the  grid.  It  is  easy  to  show  that  the  Vp  term  cannot  alter 


the  vorticity  either  since  numerically  VxVp  =  0.  Only  the  (Vp  /p.)  and 

j  J 


(l2/2  terms  can  change  vorticity,  exactly  as  dictated  by  the  physics. 

Since  the  transformation  R  is  time  reversible,  so  are  Eqs.  IV-4  -  IV-3,  so 
that  the  entire  algorithm  advances  vertex  positions  and  velocities  reversibly 
while  evolving  the  correct  vorticity  about  every  interior  vertex.  No 
numerical  generation  of  nonphysical  vorticity  can  occur,  a  rather  unique 
feature  of  both  SPLASH  and  SPLISH  among  Lagrangian  codes. 


The  pressure  p^  in  Eq.  IV-8  is  derived  from  the  condition  that  the 


new  velocities  U  .  should  be  divergence  free  at  the  new  times tep,  satisfying 


Eq.  IV-2.  The  pressure  Poisson  equation  is  derived  from  Eq.  IV-8  by  setting 


(7-  r  IJ?)^  =  0,  to  obtain  pressure  p^  such  that 


Jjt 

2 


r;  (Vp,° 


m  (v<  r  U  S  -  —  v 

i  11  j-j'i  2 


r.L?  . 


(IV-12) 


where 


Both  terms  in  Eq.  IV-12  are  simple  to  evaluate  since  the  divergence  is  taken 
over  triangle  centered  quantities.  The  paths  are  the  "surfaces"  bounding 
the  vertex  volume  of  Fig.  1,  where  the  normal  is  directed  outward  from  the 
vertex.  The  Poisson  equation  (Eq.  IV-12)  that  results  from  this  integra¬ 
tion  has  two  advantages.  First  it  is  derived  from  V 2 <+>  =  7*V<( >  as  in  the 
continuum  case.  Secondly  the  left-hand  side  results  in  the  more  familiar 
second  order  accurate  templates  (such  as  the  five-point  formula)  for  the 
Laplacians  for  homogeneous  fluids  and  regular  mesh  geometries. 

In  summary  the  finite  difference  formulas  for  SPLASH  are  derived  using 
a  control  volume  approach.  Specifications  of  pressure  at  vertices  leads 
naturally  to  the  choice  of  positioning  velocities  at  triangles.  Although 
pressure  gradients  are  constant  over  triangles,  the  resultant  algorithms  are 
expected  to  be  nearly  second  order  accurate  since  vertex  velocities  are 
derived  from  pressure  gradient  forces  through  sums  about  vertices,  which  in 
effect  centralizes  the  differences.  SPLISH  has  been  tested  extensively  on 

finite  amplitude  standing  waves  and  has  been  shown  to  be  basically  second- 

1  12 

order  accurate  by  studying  the  variation  in  period  with  mesh  size.  * 

V.  Couette  Flow 

Couette  flow  refers  to  the  circular  flow  of  a  fluid  between  two  rota¬ 
ting  coaxial  cylinders.  This  flow  is  potentially  unstable;  the  instability 
results  from  a  prevailing  adverse  gradient  of  angular  momentum. 

The  stability  of  an  ideal  fluid  in  circulatory  motion  was  first 

13 

investigated  by  Rayleigh.  Simply  stated,  Rayleigh's  criterion  says  that 
in  the  absence  of  viscosity,  the  necessary  and  sufficient  condition  for  a 
distribution  of  angular  velocity  .'Ur)  to  be  stable  is 


~~  (r2ft)2  >  0 

dr 

everywhere  in  the  interval  and  that  the  distribution  is  unstable  if  (r2&)2 
should  decrease  anywhere  in  the  interval. 

Note  that  r2ft  is  the  angular  momentum,  per  unit  mass,  of  a  fluid  element 
about  the  axis  of  rotation.  We  have  shown  in  Section  III  that  the  angular 
momentum  of  an  ideal  fluid  element  is  a  constant  of  the  motion  and  that  the 
motions  along  the  radial  and  axial  direction  may  be  treated  as  if  the 
circulatory  motion  were  absent  and  instead  a  centrifugal  force  [-pL2/2V (1/r2) ] 
were  acting  in  the  radial  direction.  Thus  we  may  associate  with  each  fluid 
element  a  "potential  energy"  pLz/2r2.  This  is  analogous  to  the  problem  of  a 
heterogeneous  fluid  in  a  field  with  a  potential  energy  proportional  to  r“2. 

The  equilibrium  is  stable  only  if  the  potential  energy  is  a  minimum;  i.e., 
the  "heavier"  fluids  are  in  regions  of  lower  potential  energy.  This  means 

that  L2  must  be  monotonically  increasing  outwards. 

3 

Taylor  extended  this  criterion  for  stability  to  account  for  viscosity, 
verified  his  calculations  experimentally  and  described  the  secondary  flow 
which  appears  after  the  onset  of  the  instability.  Viscosity  tends  to 
produce  an  angular  momentum  distribution  proportional  to  Ar2+B  for  laminar 
Couette  flow,  where  A  and  B  are  two  constants  related  to  the  angular 
velocities  of  the  inner  and  outer  cylinders.  If  this  distribution  is 
unstable,  fluid  elements  with  larger  angular  momenta  will  move  outwards 
inducing  a  secondary  flow.  Viscous  forces  will  tend  to  retard  this  motion 
but  if  a  viscous  drag  is  not  strong  enough  a  redistribution  of  angular 
momentum  will  occur.  At  the  same  time,  the  moving  solid  surfaces  will  tend 
to  re-establish  the  original  distribution  of  vorticity  and  a  steady 
secondary  flow  is  established.  This  secondary  flow  consists  of  a  regular 
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cellular  vortex  structure  in  which  closed  ring  vortices  alternating  in  sign 

are  wrapped  around  the  axis  of  symmetry. 

The  transiton  from  steady,  laminar  Couette  flow  into  fully  developed 

Taylor  vortex  flow  is  the  phenomenon  we  are  simulating  here.  To  obtain  a 

quantitative  comparison  with  theory,  we  compare  the  growth  rates  of  the 

instability  with  those  obtained  from  the  linear  theory  developed  by 
4 

Chandrasekhar . 

Following  Chandrasekhar,  we  linearize  the  equations  of  inviscid  motion 
(Eqs.  III-l  -  III-4)  by  doing  a  perturbation  expansion  about  the  stationary 
flow  solution 

u  =  eu^1)  +  e2u(2)  +  ... 

w  =  ew^1)  +  e2w(2)  +  .  .  . 

v  =  V°(r)  +  ev^1)  +  e2v(2)  +  .  .  . 

and  p  *  p°(r)  +  ep^1)  +  e2p(2)  +  .  .  . 

where  V°(r)  =  rft(r) .  Assuming  all  perturbed  quantities  vary  as  expi.i(ot-l-kz)  ] 

where  a  is  a  constant  (which  can  be  complex)  and  k  is  the  wave  number  of 
the  disturbance  in  the  z-direction,  we  have  to  first  order 


iau^1)  =  2ftv(1)  -  ~  x~  p^1^  ,  (V-l) 

p  dr 

iav^1)  +  ^  +  u^1-1  =  0  (V-2) 

iavCl)  »  -  —  p C 1 )  (V-3) 

P 

+  ikw^1^  *  0  .  (V-4) 

3r  r 


combining  Eqs.  V-l  and  V-2  and  Eqs.  V-3  and  V-4  and  eliminating  the 
pressure  between  the  resulting  two  equations  we  have  (dropping  superscripts) 


(V— 5) 


where 


[~  7~(ru) ]  *  (k2/a2)[a2  -  <t(r)]u 
dr  r  dr 


Mr)  -  f  £  (r2G)  , 


along  with  the  boundary  conditions  on  the  inner  (R  )  and  outer  (R2)  cylinders 


u(Rr)  =  u(R2)  =  0 


Since  the  numerical  code  is  Lagrangian,  we  rewrite  Eq.  V-5  in  terms  of 

4 

Lagrangian  displacement  variable  C. 
u  =  io£ 


v  =  ia?e  -  r  d?  5r 


w  =  iaC 


We  also  express  the  angular  velocity  in  terms  of  its  viscid  distribution 
il(r)  =  A  +  B/r2.  Equation  V-5  then  takes  the  form 


(DD  -k2)£  *  -  —*7-  A(A+B/r2)£ 


r  a 


(V-6) 


where 


D  =  j-  and  D .  +  - 

dr  *  dr  r 


and  the  boundary  conditions  are 


=  0  at  r  =  Rj,  R2. 


Equation  V-6  can  be  solved  explicitly  in  terms  of  Airy  functions  for 
14 

the  case  of  a  small  gap 

d  =*  (Rt  +  R2)  <<  1/2 (Rx  +  R2) . 


i 


T 


The  solutions  are  complicated  in  that  the  wavenumber  and  growth  rate  are 

linked  parametrically.  The  coupling  terms  in  the  pair  of  resulting  equations 

are  determined  from  a  characteristic  equation  expressed  as  a  ratio  of  Airy 

functions.  We  use  the  growth  rates  of  the  most  unstable  mode  as  determined 
14  4 

by  Reid  and  Chandrasekhar  to  compare  with  the  numerical  simulation. 


VI.  Numerical  Results 

The  initial  grid  for  the  Couette  flow  simulation  is  shown  in  Fig.  3. 
Ri(R2)  is  the  radius  of  the  inner  (outer)  cylinder.  zq(zn)  is  the  left 
(right)  boundary  of  the  computational  region.  The  8  coordinate  is  into  the 
page.  The  boundary  conditions  are 


u(Rx)  =  u(R2)  =  0 

and  the  system  is  periodic  in  the  z-direction 


and 


Zo  ZN+1 


w(zo)  *  W(W 


U(Z0}  *  U(W 


P<z0) 


P<? 

p(z 


N+l)  ’ 
N+l) 


The  initial  angular  velocity  is  assumed  to  have  the  viscid  distri¬ 
bution 

»l(r)  -  A  +  B/r2 

where  A  and  B  are  constants  which  depend  on  the  angular  frequency  of  the 
inner  and  outer  cylinders.  The  system  is  initially  in  equilibrium  with  a 
pressure  distribution  given  by 


p°(r,z,t=0) 


Ri 


v2(r,z) 

r 


dr  +  C 
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SPLASH 


Fig.  3  —  The  initial  computational  grid  for  the  Taylor-Couette  problem.  The 
^-direction  is  into  the  paper.  R1(R2)  is  the  radius  of  the  inner  (outer)  cylin¬ 
der.  The  r-component  of  velocity  for  vertices  1,  2  and  3  are  followed  in  time. 
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where  C  is  an  arbitrary  constant  except  in  the  case  of  a  free  surface  at 
the  inner  cylinder  when  it  is  zero. 

This  initial  pressure  distribution  is  then  perturbed  with  a  1 %  sinu¬ 
soidal  perturbation, 

p'(r,z,t=Q)  =  p°(r,z,t=0)  +  0.01p°(r,z,t=Q)  sin  kz 

where  k  is  the  wavenumber  of  the  most  unsable  mode.  The  most  unstable  mode 
is  the  mode  for  which  the  Taylor  vortices  have  a  wavelength  equal  to  twice 
the  gap  width,  i.e.,  the  vortices  are  approximately  square  in  cross-section. 

Two  important  parameters  governing  the  stability  or  instability  of 
Couette  flow  are  the  ratio  of  radii  of  the  inner  and  outer  cylinders 

n  =  R^/k^  * 

and  the  ratio  of  the  angular  frequencies  of  the  outer  and  inner  cylinders 

y  * 

The  Rayleigh  criterion  for  stability  can  then  be  written  as 

u  >  n2 

For  all  the  results  presented  here,  the  small  gap  approximation  is  valid. 

For  Case  I  we  have 

R  -  21  cm,  R  =23  cm, 

1  2 

z„  *  4  cm,  k  *  2 it /  A  ■  ir/2 
N 

and 

*  40tt  sec-1,  u  ■  1/2  , 

and  this  system  is  unstable  since  y  <  n2.  The  initial  grid  is  shown  in 
Fig.  3  where  we  have  denoted  three  vertices  which  we  will  follow  in  time. 


Figure  4  shows  the  system  shortly  after  one  full  revolution.  As  one 
would  expect  the  fluid  near  the  inner  cylinder  with  its  larger  angular 
momentum  is  being  pushed  to  a  larger  radius  and  the  fluid  near  the  outer 
cylinder  with  its  smaller  angular  momentum  is  being  pushed  to  a  smaller 
radius,  i.e.,  the  "heavy"  fluid  is  falling  and  the  "light"  fluid  is  rising. 
Note  the  large  number  of  vertices  that  have  been  added  near  the  boundaries 
to  preserve  the  resolution  there. 

The  time  evolution  of  the  r-component  of  velocity  for  the  three 
aforementioned  vertices  is  shown  in  Fig.  5.  The  growth  rates  for  the  three 
vertices  are  equal  for  nearly  a  full  revolution  at  which  point  vertex  1 
slows  down  and  vertices  2  and  3  speed  up.  The  growth  rate  in  the  linear 
regime  yc  =  213.09  s-i  is  in  very  good  agreement  with  the  predicted  growth 
rate  y^  =  215.26  s-1  obtained  from  linear  theory. 

Vertices  2  and  3  are  speeding  up  because  they  are  becoming  entrained 
between  two  very  large  counter-rotating  vortices.  Vertex  1  is  nearing  the 
inner  cylinder  and  being  deflected  in  the  z-direction.  This  is  shown 
clearly  in  Fig.  6  where  we  have  plotted  the  pathlines  of  the  flow.  The 
plus  signs  are  the  most  recent  positions  of  the  vertices  and  the  dots  are 
the  positions  at  three  previous  equal  periods  of  time.  The  plus  signs  can 
be  regarded  as  arrow  heads.  We  see  the  development  of  two  large  counter¬ 
rotating  Taylor  vortex  cells.  The  wavelengths  of  the  vortex  cells  is  4  cm 
as  predicted  by  the  theory. 

Figures  3,  4  and  6  can  be  directly  compared  to  the  experimental  work 
of  Donnelly  and  Fultz^  (see  Fig.  7).  They  ejected  dye  from  the  inner 
cylinder  and  took  a  remarkable  series  of  photographs  showing  the  transition 

from  laminar  flow  to  fully  developed  Taylor  vortex  flow.  Their  work  is  also 

4 

illustrated  in  Chandrasekhar. 
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PATHL  I  NES 


rf.lJMlT 


Fig.  6  —  Pathlines  for  Case  I.  The  +  signs  are  the  most  recent  positions  of  the 
vertices  and  the  dots  are  the  positions  at  three  previous  equal  periods  of  time. 
The  center  vortex  has  a  clockwise  rotation.  Vertices  1,  2  and  3  are  noted. 


The  onset  of  instability  with  the  outer  cylinder  at  rest,  n  =  0,  Pc  =  4-4'Jl  soo: 

(a)  laminar  flow,  P  ==  4-500  soc  ;  (b)  beginning  of  radial  motion  at  P  —  4-483  sec ;  (c)  Appear¬ 
ance  of  cells  with  the  outer  cylinder  at  rest;  cells  at  marginal  stability,  P  —  Pe  =  4-466  soc; 

(d)  Appearance  of  cells  with  cylinders  rotating  in  the  samo  direction:  P  —  Pe  =  3-844  sec, 

p  =  0-1164. 

Fig.  7  —  Photographs  showing  the  transition  from  laminar  Couette  flow  to  fully  developed 
Taylor-vortex  flow  with  the  outer  cylinder  at  rest.  Dye  is  ejected  from  the  inner  cylinder. 
It.J.  Donnelly  and  D.  Fultz,  Proc.  Roy  Soc.  A  258,  101  (1960).  Used  by  permission. 
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For  Case  II,  the  outer  cylinder  is  at  rest,  y  =  0,  and  the  system 
length  has  been  increased  to  12  cm.  R^ ,  R^ ,  and  k  are  the  same  as  for 
Case  I.  Since  the  wavelength  of  the  perturbation  is  4  cm,  six  Taylor  cells 
should  develop.  This  is  illustrated  in  Fig.  8  where  we  have  plotted  contours 
of  constant  circulation.  The  plus  sign  indicates  flow  in  the  clockwise 
direction. 

The  case  for  which  the  cylinders  are  rotating  in  the  opposite  direction 
is  particularly  interesting.  Theory  predicts  that  only  the  inner  region  of 
fluid  should  be  unstable  as  it  is  only  in  this  region  that  the  angular 
momentum  is  decreasing.  For  Case  III  we  have 
R1  =  41  cm,  R2  =  45  cm, 
z^  =  8  cm,  k  =  it 

=  40 ns-1,  y  =  -0.87. 

This  choice  of  y  gives  zero  angular  momentum  at  the  center  of  the  gap.  The 
constant  circulation  contours  are  shown  in  Fig.  9.  As  predicted,  four 
counter-rotating  vortex  cells  appear  in  the  inner  half  of  the  fluid  while 
the  outer  half  of  the  fluid  remains  stable.  This  figure  can  also  be 
directly  compared  with  the  photographo  of  the  experimental  work  of  Donnelly 
and  Fultz1'*  (Fig.  10). 

The  numerical  results  for  all  the  cases  are  summarized  in  table  form 
in  Fig.  11.  The  computational  growth  rates  are  in  very  good  agreement  with 
the  theoretical  growth  rates  with  errors  on  the  order  of  1%.  Note  that  for 
u  =  1  the  fluid  is  stable  and  y  is  the  oscillation  frequency  of  the 
perturbed  velocity. 
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NUMERICAL  RESULTS 
COUETTE  FLOW 
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Fig.  11  —  Comparison  of  computational  growth  rates  with  growth  rates  obtained 
from  linear  theory  for  various  initial  conditions 
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VII.  Summary  and  Conclusions 

We  have  developed  a  2  1/2  dimensional,  Lagrangian,  hydrodynamic  model 
designed  for  Che  simulation  of  transient  rotational  flow  phenomena.  The 
model  uses  as  a  finite  difference  mesh  a  general  connectivity  triangular 
grid.  The  advantages  of  this  model  are  numerous: 

1)  Complicated  geometries,  interfaces  and  free  surfaces  can  be 
treated  with  a  minimum  of  difficulty. 

2)  The  resolution  across  the  mesh  is  highly  variable. 

3)  The  mesh  can  be  restructured  to  preserve  the  numerical  accuracy 
of  the  simulation. 

4)  No  numerical  vorticity  is  generated. 

5)  The  gradient  operator,  as  a  triangle  function,  is  exact  in  the 
linear  approximation. 

6)  The  integral  operator,  as  a  vertex  function,  is  exact  in  the  linear 
approximation. 

We  have  applied  this  model  to  the  study  of  the  transition  of  laminar 
Couette  flow  to  Taylor  vortex  flow  with  a  high  degree  of  success.  The 
computational  growth  rates  are  in  excellent  agreement  with  the  theory. 

For  the  results  presented  here  the  code  has  not  been  run  long  enough 
to  achieve  steady-state  Taylor-vor tex  flow.  However,  with  respect  to 
the  transition  to  Taylor-vortex  flow  we  find  that  a  vortex  signature  appears 
very  early  in  the  run  when  the  system  is  still  in  the  linear  regime.  By 
this  we  mean  that  contours  of  constant  circulation  show  uniformly  spaced 
vortices  when  the  perturbed  velocities  are  on  the  order  of  10~6  cm/sec. 

[V.  0(103  cm/sec)].  These  vortices  then  increase  in  strength  but  maintain 


their  shape  and  spacing. 


In  these  calculations  we  have  perturbed  the  system  at  only  one  wave¬ 
length  corresponding  to  the  predicted  wavenumber  for  the  steady-state  Taylor 
vortices.  Although  the  wavenumber  is  therefore  not  expected  to  change,  it 
is  surprising  to  find  that  the  flow  throughout  the  cylinder  gap  is  estab¬ 
lished  during  the  linear  regime  in  exactly  the  nonlinear  flow  pattern.  As 
evidenced  by  the  case  of  the  counter-rotating  cylinders,  this  flow  is  not 
evident  a  priori.  The  transition  to  Taylor-vortex  flow  from  laminar 
Couette  flow  is  strickly  speaking  nonexistent —  the  linear  flow  is  only 
strengthened.  Whether  any  consequences  of  this  simple  transitioning  are 
evident  in  the  more  complicated  cases  of  viscid  flow  a  d  the  perturbations  of 
many  wavenumbers  will  be  investigated  in  future  calculations. 

Two  code  modifications  are  necessary  for  these  calculations  to  be 
compared  with  experiments;  the  addition  of  viscosity  (see  Appendix  A)  and 
the  improvement  of  some  grid  restructuring  algorithms.  In  order  to 
preserve  numerical  accuracy  when  the  grid  becomes  highly  distorted  vertex 
additions  and  deletions  must  be  made.  We  have  now  developed  vertex 
addition/deletion  algorithms  which  conserve  divergence,  curl,  linear 
momentum,  angular  momentum  and  vorticity  generation.  These  algorithms  are 
now  being  incorporated  into  the  code  to  allow  the  calculation  to  proceed 
further  into  the  fully  nonlinear  regime  while  preserving  second-order 
accuracy.  These  vertex  addition/delection  algorithms  do  not  affect  the 
results  presented  here. 

The  inclusion  of  viscosity  and  the  new  restructuring  algorithms  will 
enable  a  direct  comparison  between  the  computational  results  and  the 
experimental  results  to  be  made.  This  would  entail  calculating  the  torques 
required  to  maintain  the  cylinders  in  motion  and  the  critical  Taylor  number, 
the  ratio  of  the  destabilizing  centrifugal  force  and  the  opposed  radial 
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Appendix  A 


Consider  the  viscid  equations  of  motion  written  in  our  pseudo- 
Cartesian  format 

4“  .  -  I  Xl  u  f2J./3u)+  J_  /3u  +  Sw)  +  2/3u  _  u\]  ( 

dt  p  ar  +  r  +o[23rl3r/  3z  Uz  +  3r/  +  r(3r  r/J*  (A  1} 


and 


dw 

dt 


dL 

dt 


=  _Ii2  +  M./2—  [— 1  +  -  —  Fr/—  +  — \11  (A-2) 

p  32  p  1  32  1.3 zj  r  3r  Lr\3z  3r/Jj  '  ' 


f 32v  ,  3  /  3v  v\  2  /  3v  v\1 

rpL3z2  3r  \  3r  -  r )  r  \  3r  r/J 


The  coefficient  of  viscosity  p  is  considered  a  constant. 

A-2  can  be  combined  to  give 

du  - 

5J---7P-  1/21^ -£±7x0*)  , 

where  uj_  is  the  theta  component  of  vorticity 


(A-3) 

Equations  A-l  and 

(A-4) 


uj  =  u)ae  =  Vxu 
o  y 

u)  is  a  vertex  function  which  is  easily  determined  by  calculating  the  circula- 

0 

tion  about  a  vertex 


u»,  =  ~  §  u’dl 

8i  Ai  “  ~ 


Equation  A-3  can  be  cast  into  the  form  of  a  diffusion  equation  for  the 
angular  momentum. 

^-ryV-^VL)  .  (A-5) 

The  finite  differencing  of  this  equation  is  discussed  in  detail  in  Ref.  10. 
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